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Abstract 

We study the distribution of first passage time (FPT) in Levy type of anoma- 
lous diffusion. Using recently formulated fractional Fokker-Planck equation 
we obtain three results. (1) We derive an explicit expression for the FPT dis- 
tribution in terms of Fox or H- functions when the diffusion has zero drift. (2) 
For the nonzero drift case we obtain an analytical expression for the Laplace 
transform of the FPT distribution. (3) We express the FPT distribution in 
terms of a power series for the case of two absorbing barriers. The known 
results for ordinary diffusion (Brownian motion) are obtained as special cases 
of our more general results. 
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I. INTRODUCTION 



For a stochastic process, the first passage time (FPT) is defined as the time T when the 
process, starting from a given point, reaches a predetermined level for the first time, and is 
a random variable [0-^. In this paper, we consider the distribution of FPT in the context 
of anomalous diffusion. By anomalous diffusion we mean a process where the mean square 
displacement of the diffusive variable X{t) scales with time as ~ P with < 7 < 2. 

If 7 = 1 we get ordinary diffusion. There are different mechanisms for generating anomalous 
diffusion. Our focus will be on the continuous time random walk (CTRW) where the 
waiting time is a Levy type variable obeying certain power law distribution. We subsequently 
refer to the process as the Levy type of anomalous diffusion. Recent work shows that in 
the generalized diffusive limit the probability density function of the CTRW is described 
by a fractional Fokker-Planck equation (FFPE) Using this equation, together with 

appropriate initial and boundary value condition, we derive the exact solution for the FPT 
density function in terms of Fox or H functions P,pi]| when the anomalous diffusion has zero 
drift. For the nonzero drift case we derive an expression for the Laplace transform of the 
FPT density function. This Laplace transform is then used to obtain the mean and variance 
of the FPT. Finally, we derive an explicit expression for the FPT density function for Levy 
type anomalous diffusion with absorbing barriers at finite distances from the origin. 



II. LEVY TYPE ANOMALOUS DIFFUSION 



For an ordinary diffusive process, the mean square displacement scales with t as {X^(t)) ~ 
t in the large t limit. If (X^(t)) ~ P, we have anomalous diffusion, with < 7 < 1 and 
1 < 7 < 2 corresponding to subdiffusive and superdiffusive cases, respectively. Anomalous 
diffusion with drift X^{t) is defined by 

X^it)= fit + X{t). (2.1) 

For further discussions and examples of anomalous diffusions in physical, chemical and 
biological systems see PJ^. For completeness we give in what follows a brief introduction 
to CTRW and the limiting process leading to the establishment of the FFPE. 



A. Details of the Process 

We start with a one dimensional continuous time random walk described by the following 
Langevin equation: 

rlX °° 

— = ^F,5(t-t,). (2.2) 

ai .^-^ 

Here the random walker starts at x = at time to = 0. Subsequently, the random walker 
waits at a given location Xi for time ti — before taking a jump Yi which could depend on 
the waiting time. The waiting time m > and the jump size y (—00 < y < 00) are drawn 
from the joint probability density function (j){y,u). The waiting time distribution ip{u) is 
given by 
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i,{u)= / dyct){y,u). (2.3) 

J —oo 

The process is non Markovian if '4){u) is a non-exponential distribution since the probabihty 
for the next jump to occur depends on how long the the random walker has been waiting 
since the previous jump. But CTRW is non Markovian in a special way since it does not 
depend on the history of the process prior to the previous jump. 

It can be shown |^ that the probability distribution W{x,t) for the CTRW is related to 
(f){y,u) and '?/'(«). This relation is simpler to write down in the Fourier-Laplace transform 
space. Denoting the Fourier-Laplace transform of W{x, t) and (f){y, u) by W{k, s) and 0(/c, s) 
respectively (where k is the Fourier transform of the space variable and s the Laplace 
transform of the time variable) we have 

W(Ks)^\^^^. (2.4) 

Here ip{s) is the Laplace transform of ip{u). 

It can be further shown that, depending on the specific form of (f){y,u), the CTRW can 
produce both subdiffusive (0 < 7 < 1) and superdiffusive processes (1 < 7 < 2) as well as 
ordinary diffusion (7 = 1) P,p!T|. For example, consider 

0(j/, u) = -jL= exp[-yV2(T'] ~ ^'^J'' , (2.5) 

where y and u are decoupled with y being a Gaussian variable with zero mean. For 1 < a < 2, 
the corresponding CTRW is characterized by a subdiffusive process with 7 = a — 1, and for 
a > 2, one gets ordinary diffusion with 7 = 1. This is demonstrated as follows. 
The waiting time distribution ip{u) is given by [cf. Eqs. ( |2.5{ ) and ([^.3D] 

(1 + ?i/r)" 

Even though this is not a Levy stable distribution, it belongs to the domain of attraction 
1^ of a one-sided stable Levy law. It has infinite mean for 1 < a < 2, the range we are 



interested in (it has a finite mean for a > 2). We call u a "Levy type of variable" for want 
of a better name. The Laplace transform 'il){s) of ip{u) is ||T3| 

^{s) = (a - l)(rs)"-^r(l - a, rs)e"^ (2.7) 

The Fourier-Laplace transform (j){k, s) of (j){y, u) in Eq. ( p. 51 ) is given by [|13| 

0(A;, s) = exp(-a^A;V2)^(s). (2.8) 



To show that the CTRW characterized by Eq. ( p.5| ) exhibits anomalous diffusion, we 
need to demonstrate that {X'^(t)) ~ for large t (0 < 7 < 1). By Tauberian theorem |]14 
this is equivalent to (X^(s)) ~ l/s^+^ for small s (0 < 7 < 1). But we have the result 



2 



k=0- 



(2.9) 
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Thus we need the expression for W{k, s) in the hmit s — > 0. Since W{k, s) is a function of 
ip{s) and 4>{k,s) [cf. Eq. (|2.4|) ], we first consider ipi^s). We have 

r(l - a, rs) = r(l - a) - E ^ J , . • (2.10) 



n=0 



n 



!(1 - a + n) 



As s 0, 



r(2-«) (rs)i-" (rs)2-" 

r 1 - «, ^ + + 2.11 

a — 1 a — I 2 — a 



Therefore ip{s) as s — > is given by 



^p{s) ^l-T{2-a){Tsf-\ l<a<2, (2.12) 
^ 1- (2a-3)rs/(a-2), a > 2. (2.13) 



Substituting this in the Eq. ( p.4|) we have [cf. Eq. (|2.8|) ] 

sl-[l- r(2 - a)(rs)"-i] exp(-a2F/2)' ^ < « < ^' ^^'^"^^ 

^^'/(^ - 2) (2 1^) 

^ sl-[l-aTs/{a-2)]exp{-a^ky2y ^ ' 



Using this in Eq. (|2.9|) and evaluating the second derivative at A; = we obtain in the hmit 



<^(-)^- r(2-.).--^ ;^' 1<«<2, (2.16) 

a>2. (2.17) 
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Thus we have shown that the CTRW characterized by Eq. ( p.5| ) is a subdiffusive process 
for 1 < a < 2 with 7 = a — 1, and for a > 2, one gets ordinary diffusion with 7 = 1. For 
a = 2, by taking proper hmits, one can show that we again obtain ordinary diffusion. 

We verify the above resuh numerically by simulating the CTRW process. For the sake 
of numerical efficiency, we replace the waiting time distribution tp{u) in Eq. (p^.6D by the 



Pareto distribution [IG 



ij{u) = 0, u<T, (2.18) 
(a-l)r°-i 

-, u>T, (2.19) 



This approximation is well justified for small values of r. Numerically obtained mean square 
displacement is compared with the theoretical prediction in Figure 1 for a = 1.5 (i.e. 7 = 
0.5), a = 1.0, r = 10"^ and = 3.5 x 10~^. We note that the numerical simulation is in 
excellent agreement with the theoretical prediction. 

If, on the other hand, (f){y, u) is given by the coupled distribution p 
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^( ^ I 11/ A (/g-l)/^ 

u) = -5{u/t - \y\/a)j^^^-j-^^, 



(2.20) 



where 2 < /? < 3 and 5{-) is the Dirac delta function, the CTRW describes a superdiffusive 
process with 7 = /? — 1. This can be seen as follows. 

Proceeding as before, the waiting time distribution ip{u) is given by [cf. Eqs. ( 2.201 ) and 

:3Di 



ip{u) 



(/3-l)/r 



(2.21) 



It has a finite mean but infinite variance for 2 < /3 < 3. Its Laplace transform in the limit 
s — > is given by 



V^(s)^l-(2/5-3)rs/(/3-2) 



(2.22) 



The Fourier-Laplace transform 0(/c, s) of 0(?/, u) in Eq. ( 2.20| ) is given by (in the limit s — 
and -C s) 



(2.23) 



Substituting ipi^s) and 0(fc, s) in Eq. (|2.4| ) we get 
W(fc,s) 



s + 



(/3-l)(/3-2)r(3-/3) ,,^ ,_3 
2(2/3-3) 



(2.24) 



Using this in Eq. 
s ^ 



and evaluating the second derivative at = we obtain in the limit 

(/?-i)(/3-2)r(3-/3) 1 



{X\s)) 



(2/5 - 3) 



,5-/3- 



(2.25) 



By Tauberian theorem |T^, this is equivalent to (X^(t)) ~ Thus the CTRW charac- 

terized by Eq. ( |2.20| ) is a superdiffusive process for 2 < /3 < 3 with 7 = 4- /3 (1<7<2). 

Note that there are several other possible forms for 0(|/, u) which also give rise to sub- 
and superdiffusive behavior. See for a detailed discussion. 



B. Fractional Fokker-Planck Equation for Levy type anomalous diffusion 

In the previous subsection, we had described CTRW processes that give rise to Levy type 
anomalous diffusion. However, it is difficult to derive analytical results directly from the 
process. It is more convenient to work in the general framework of Fokker-Planck equations 
One can go from the CTRW process to a fractional Fokker-Planck equation (FFPE) 

by taking the generalized diffusion limit. This limit is analogous to the regular diffusion 
limit that one considers to derive the regular Fokker-Planck equation from a random walk 
with (X^(t)) ~ t for large t. In the regular diffusion limit, one lets cr, r ^ such that a^/r 
is maintained a constant. For a CTRW where (X^(t)) ~ t'^ (0 < 7 < 2), the generalized 
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diffusion limit is obtained by taking the limit o", r — such that cf'^/t'^ is maintained a 
constant. 

Thus to obtain a FFPE from the CTRW process, we need to take the limit a, r ^ 0. 
We will take this limit for the Fourier-Laplace transform W{k, s) of W{x, t) and then invert 
the transform to obtain the FFPE. First consider the subdiffusive CTRW characterized by 
Eq. (|2.5|) . In the previous section, we have already derived an expression for W{k, s) in the 
limit s ^ (which is equivalent to the limit r that we now require): 

^) " n rn^^\/\y ( 2.2 „y (2-26) 
s 1 — [1 — r(l — 7)(rs)^J exp(— cr"'A;"'/2) 

Here we have used 7 = 0; — 1 (0<7< 1) instead of a since it is the physically relevant 
quantity. Now we take the further limit a ^ such that a"^ /2T{1 — 7)r'^ = 7^ is a constant. 
K is called the generalized diffusion constant. We then obtain 

This can be rewritten as 

W{k, s)-\ = -Kk'^s-m{k, s). (2.28) 

To take the inverse Fourier-Laplace transform of the above equation, we need the inverse 
Laplace transform of s~'^W{k, s). This is given by the Riemann-Liouville fractional integral 
QDi^W{k,t) which is defined as 0JT|] 



oD^-^Wik, t) = -^ f dt' (t - ty-'Wik, t'), 7 > 0. (2.29) 

r(7) Jo 

This result enables us to take the inverse Fourier-Laplace transform of Eq. ( |2.28|) giving 

Wix,t)-Wix,0)= K oD;^—Wix,t), 0<7<1. (2.30) 

Here we have incorporated the initial condition W{x, 0) = 6{x). 

Next we consider the superdiffusive CTRW process characterized by Eq. ( [^.2(J| ). We 
have already derived an expression for W{k, s) in the limit fc, s — > (which is equivalent to 
the limit a, r — that we now require): 



W{k,s) 



(3-,)(2-,)rh-i) 

2(5-27) ^ ' 



(2.31) 



This expression is valid in the regime a ^ r. Taking the generalized diffusion limit we again 
get Eq. (|;28D but now K = (3-7)(2-7)r(7-l)aV2(5-27)r^ and 1 < 7 < 2. Therefore, 
this superdiffusive CTRW process also gives Eq. ( p.30|) with the above K and 1 < 7 < 2. 

To summarize, a class of Levy type anomalous diffusion can be described by the following 
FFPE in the generalized diffusion limit: 
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rj2 

W{x,t)-W{x,0)= K oD;^—W{x,t), 0<7<2. (2.32) 

Here K has different expressions for < 7 < 1 and 1 < 7 < 2 as described above. We 
note that there are other CTRW processes which may give rise to a different FFPE. In this 
paper, we consider only the FFPE given in Eq. (|2.32|) . The above FFPE describes Levy 
type anomalous diffusion with zero drift. For anomalous diffusion with drift /i, the above 
FFPE can be generalized to give [^]: 

Wix, t) - W{x, 0)= K ,D;^—Wix, t)-ii ,D;^—W{x, t), (2.33) 

where < 7 < 2. We also note that for 7=1, the above FFPE reduces to the regular 
Fokker-Planck equation describing Brownian motion. 

We now show that this FFPE gives the correct moments. The natural boundary condi- 
tions for the FFPE are: 

W{-oo, t) = W{oo, t) = 0, W{x, 0) = 6{x). (2.34) 

First we compute (X^). Multiplying the FFPE given in Eq. (|2.33| ) by x and integrating 
from X = —00 to X = 00 we get 

{X,) = oD;' r dx W{x, t)-K oA"" r dx ^^pll^ (2.35) 

J —oc J —oo ox 

where we have performed integration by parts and assumed that xW{x, t) = and x ^^g^'^^ = 
at a; = ±00. Evaluating the remaining integrals using the natural boundary conditions 
and the normalization integral J^^dxW{x,t) = 1, we get 

{X^) = fi oD^\l) = fit. (2.36) 

This is the expected result. 

Next we compute Var[X^]. Multiplying the FFPE in Eq. (|2.33|) by x'^ and integrating 
we get 

(Xj) = 2/i oA"^ r dxxW{x,t) - 2K oA'^ H dxx ^^^^'^\ (2.37) 
J —00 J— 00 ox 

where we have again used integration by parts and assumed that x'^W{x,t) = and 
x^^^^^^ = at X = ±00. Evaluating the remaining integrals using the fact that 
rfa;a;iy(a;, t) = (X^) = fit and the boundary and normalization conditions we obtain 

{XD = 2/i2 ,D-\t) + 2K oA""(l)- (2.38) 
But oDi\t) = ty2 and oA"^(l) = tyT{-f + 1) 0. Therefore 

The variance Var[X^] is given by (X^) — (X^)^. Hence 

Thus we see that our FFPE gives the correct moments. 
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III. FIRST PASSAGE TIME PROBLEM FOR LEVY TYPE ANOMALOUS 

DIFFUSION 



We now formulate the first passage time problem for the FFPE given in Eq. (p.33| ) 
describing a Levy type anomalous diffusion with drift. Consider a stochastic process X{t) 



with X(0) = 0. The first passage time (FPT) T to the point X = a is defined as |T9 



T = inf{t : X{t) = a}. (3.1) 

We would like to obtain the probability density function for T. This is the first passage time 
problem. 

For a process described by Fokker-Planck equations, the problem of obtaining the FPT 
density function can be recast as a boundary value problem with absorbing boundaries [|I|. 
In our case, to obtain the FPT density function, we first need to solve Eq. ( p.33|) with 



absorbing boundaries at x = — oo and x = a, where a is the predetermined level of crossing, 
with the initial condition W{x,0) = 6{x) An equivalent formulation, due to symmetry, 
is to solve Eq. ( [^.33D with the following boundary and initial conditions: 

W{0,t) = 0, W{oo,t) = 0, W{x,0) = 5{x~ a), (3.2) 

where x = a is the new starting point of the CTRW, containing the initial concentration of 
the distribution. The equivalence is easily seen by making the change of variables x a — x 
in Eq. ( |2.33| ). The only change is in the interpretation of fi. Now /i < corresponds 



to drift towards the barrier. This latter formulation makes the subsequent derivation less 
cumbersome. Once we solve for W{x,t), the first passage time density f{t) is given by 

f{t) = -TJ dxW{x,t). (3.3) 
at Jo 

In the following subsections, we obtain the FPT density function for Levy type anomalous 
diffusion under different conditions. 



A. Zero Drift Case 

In this section, we obtain an explicit solution for the FPT density for a zero drift Levy 
type anomalous diffusion in terms of H-functions P, p!0| , |20[| . Asymptotic expressions for the 



FPT density in various limits are obtained. Numerical simulations are performed to confirm 
the above results. 

Setting /i = in Eq. ( |2.33| ) we obtain 

W{x, t) - W{x, 0)= K ,D;^—W{x, t). (3.4) 

Taking into account of the boundary and initial conditions [cf. Eq. (|3.2| )] we are led to the 
following expansion for W{x,t) [pl|] 



2 f°° 

W(x,t) = — dk sin kx sin ka A(k,t), (3.5) 

TT JO 
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with A{k,0) = 1. To determine the unknown function A{k,t), we substitute the above 
expansion for W{x,t) in Eq. ( p.4|) and, after straightforward algebra, obtain A{k,t) — 1 = 
—Kk"^ QD^"'A{k,t). Taking the Laplace transform with respect to t, we have 



A{k,s) 



s + k^Ks^-"/' 



(3.6) 



where A{k,s) is the Laplace transform of A{k,t). Here we have applied the result \TE\ that 
the Laplace transform of qD^^ A{k^t) is A{k, s) / s'^ . 
Inverse Laplace transform of Eq. ( |3.6| ) yields |^| 



A{k,t) = E^{-k^Kf<), 



(3.7) 



where E^(z) is the Mittag-Leffler function |2^. The Mittag-Leffler function can be defined 
by the following power series expansion: 



n=Q 



r(n7 + 1) 



(3.8) 



When 7 = 1, the Mittag-Leffler function reduces to the usual exponential function e^. 
Substituting Eq. into Eq. (|]|) we get 



(3.9) 



2 

W{x,t) = — dksmkx sin ka E-^{-k'^Kf). 

TX Jo 



To proceed further, we introduce the Fox or H-function P, p!0| , pO| which has the following 
alternating power series expansion: 



Tjm,n 



(ttj, Aj)j=i^ 



E E 

1=1 k=0 



k\B, 



UT=l,m r(^J - BjSik) Ur=l r(l -ar + ArSik) 
nLm+1 r(l -bu + BuSik) nLn+1 ^{(^v - Ksik) 



(3.10) 



where sik = ipi + k)/Bi and an empty product is interpreted as unity. Further, m,n,p,q 
are nonnegative integers such that < n < p, 1 < m < q; Aj, Bj are positive numbers; 
aj, bj can be complex numbers. We will denote the H-function by H^^^{z) for the sake 
of notational simplicity wherever this does not lead to any confusion. The H-function has 
several remarkable properties, some of which are listed in the Appendix. 

Returning to our problem, by comparing the series expansion [cf. Eq. ( |3.8| )] of the 
Mittag-Leffler function E^{z) with that of the H-function [cf. Eq. ( |3.10|) ], we see that 



E^ 



bH 



(0,1) 

(0,1), (0,7) 



(3.11) 



Substituting this in Eq. 



W{x,t) = - 

IT Jo 



9f) we obtain 
dksmkx sm ka I k'^Kt"' 



(0,1) 

(0,1), (0,7) 



(3.12) 
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Letting k' = k^KVy/"^ the above equation becomes 



W{x,t) 



dk'sink'x sink'a H^'^ (^{kT 



(0,1) 

(0,1), (0,7) 



7r(irn)i/2 

Using Property 5 [Eq. ( |A4D ] of the H-functions to replace [k'Y by k' we get 



(3.13) 



W{x,t) 



1 11 f 

— r-jT / dk' sin k'x sin k'a '2 k' 

tf)^/'^ Jo ' \ 



(0, 1/2) 

(0,1/2), (0,7/2) 



(3.14) 



The above equation can be rewritten as follows making use of the standard trigonometric 
identity 2 sin k'x sin k'a = cos k'{x — a) — cos k'{x + a): 



W{x,t) 



1 f°° 11(1 
^^^yl2 J dk' [cos k' {x — a) — COS k'{x + a)] H 12 \k' 



(0, 1/2) 

(0,1/2), (0,7/2) 



(3.15) 



The above Fourier cosine transforms can be solved by successive applications of a Laplace 
and an inverse Laplace transform (a technique pioneered by Fox p3[ for solving a wide 
variety of integral transforms) to give pO[ 



W{x,t) 



\x — a 



2\x — a\ 
1 

2{x + a 



3,3 



2,1 f X + a 

3,3 



-H. 



(1,1/2), (1,7/2), (l,l/2)\ 
(1,1), (1,1/2), (1, 1/2)^1 

(1,1/2), (1,7/2), (l,l/2)\ 
(1,1), (1,1/2), (1,1/2) ]• 



(3.16) 



Now, applying Property 2 [cf. Eq. (|Al|)] of the H-functions to reduce the order of our 
H-function, we get 



W{x,t) 



1 



7H. 



2,0 
2,2 



\x — a 



2.0 / X -\- a 



-H. 



2{x + a) '•'{{Kt^y/^ 



(1,7/2), (l,l/2)\ 
(1,1), (1,1/2)^ 
(1,7/2), (l,l/2)\ 
(1,1), (1,1/2) j 



(3.17) 



Applying Properties 1 and 3 [cf. Eq. ([X^)] of the H-functions, we can further reduce the 
order of our H-function: 



Wix,t) = -^Hl' 



1 


X — a 




(l,7/2)\ 
(1,1) ) 


\{Ktyy/^ 



1 rri.o ( 3: + a 

'2{x + a) ^^^[{Kt^y/^ 



(l,7/2)\ 
(1,1) j 



(3.18) 



Applying Property 6 [cf. Eq. ( |A5D ] of the H-functions with z = \x — a\ (or z = x + a) and 
p = — 1 we finally obtain 



W{x,t) = 
1 

2{Kriy/^ 



-"1,1 



\x — a\ 



{Kt^y/^ 



(1-7/2,7/2) 

(0,1) 



1,0 



X + a 

{Kt^y/^ 



(1-7/2,7/2) 

(0,1) 



(3.19) 
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Substituting Eq. ( p.l9| ) into Eq. ( p.3|) we have 



fit) = - 



d 
di 



1 



1,0 



d 

^dT 



2(irt7)i/2 / 



dx 



X — a 



J 7^1,0 a; + a 
ax ill 1 



(1 -7/2, 7/2) y 
(0,1) )_ 
1 -7/2, 7/2) \ 
(0,1) ' 



Defining z = {x - a)/{KT^y/^, z' 

d 1"°° 

m = -4 



X + a)/{KT^y/^, we obtain 



(3.20) 



dz I \z 



dt J-a/{Kt-iy/^ 

d /■°° 



dz' Hl;^^ { z' 



dt Ja/iKtf)^/^ 



(1-7/2,7/2) 
(0,1) 

(1 -7/2, 7/2) \ 
(0,1) j 



(3.21) 



Evaluating the above equation (which is easily done since the only t dependence is in the 
limits of the integrals), we obtain the following expression for the first passage time density 



fit) 



07 



2XV2t(2+7)/2 



-"1,1 



(irt^)V2 



;i-7/2,7/2)\ 

(0,1) )■ 



(3.22) 



We mention that this result has appeared earlier in a short communication . It should be 
noted that H-functions were first used in the context of probability distributions by Schneider 



25|. They have also been used to express solutions of fractional diffusion equations p6 



The series expansion of the H-function in Eq. (|3.22|) [cf. Eq. (|3.1CI|) ] is 

07 ^ {-a/{Kr'y/y 



f(t) = '11 V 

■I ^ ^ 2KyHi^+^)/^ A;!r(l - 7/2 - A;7/2) ' 



(3.23) 



This turns out to be also the series expansion of the Maitland's generalized hypergeometric 
function or the Wright function o'ipi [H]- Thus, an alternative expression for f{t) is 



fit) 



07 



2i^l/2^(2+7)/2 



oi'i 



(1-7/2,-7/2)' (Kt7)l/2 



(3.24) 



It is difficult to work with the series expansion of f{t) given in Eq. ( |3.23|) due to the 
presence of the gamma function with large negative arguments. We get around this as 
follows. From the properties of gamma function ||15| we have 



r(l - 7/2 - A:7/2) = J±±^Ti-ik + 1)7/2). 



(3.25) 



Using the following relation [15 



Ti-z) 



-TT 



sin(7r2;)r(z + 1) ' 



(3.26) 



Eq. (|3.25|) can be rewritten as 
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r(l-7/2-A;7/2) 



7i{k + 1)7 



1 



2 smp+ l)7r7/2]r[l + (A; + 1)7/2] 
n 



sinp + l)7r7/2]r[(fc + 1)7/2]' 
Substituting this in Eq. ( p.23|) , we obtain 



(3.27) 



fit) 



07 



2vriri/2t(2+7)/2 



E 



'-a/{Kt 



sin[(A; + l)7r7/2]r[(A; + 1)7/2] 



(3.28) 



Note that for regular Brownian motion (7 = 1), the above expression for f{t) reduces to 

[-a/^Kif sin[(A; + l)7r/2]r[(fc + l)/2] 



fit) 



fc=0 



But [15 



T{k + l) 



T{k + l] 



4^/'r[(fc + i)/2]r(i + fc/2) 



(3.29) 



(3.30) 



Further, 



sin[(A; + l)7r/2] = 0, odd, 

= i-lf'^ k even. 



Substituting these relations in Eq. ( |3.29| ) and letting n = A;/2 we get 

-a/y/Kif''V[{2n + l)/2](-l)" 



fit) 



E 



n=0 



Simplifying this we finally obtain 

fit) = 



(4)"r[(2n+ l)/2r(n + 1) 



E 



(3.31) 



(3.32) 



nl 



■.exp[-a'^/4:Kt]. 



(3.33) 



This is the expected inverse Gaussian distribution for the FPT density function of the 



ordinary Brownian motion with zero drift [|19 



Next we consider the asymptotic behavior of the FPT distribution for large values of t. 



Refer to Eq. (|3]2|). Let z = a/iKVf^. It is known that for small z, Hlf 

|^|6i/-Bi _ since 61 = and Bi = 1. Therefore, the FPT distribution fit), for large t, is 
characterized by the power law relation 

fit) ~ (3.34) 

which becomes the well known —3/2 scaling law for the ordinary Brownian motion. We 
make two comments. First, the above power law behavior has been observed earlier by 
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Balakrishnan p8|] for subdiffusive processes (0 < 7 < 1) using a different method. Using 
our method the same scahng law is shown to be apphcable also to superdiffusive processes. 
Second, from Eq. ( p.34| ), we see that the mean first passage time and all higher moments of 
the FPT distribution are undefined for < 7 < 2. 

Next we consider the asymptotic behavior of f{t) for small t (i.e. large z where z = 
a/^KVY/"^). It is known p|j2l that for large z, 



where 



exp 



(3.35) 



« = E - E «i + - ^ + i)/2, 

j=i 1=1 

q p 

/^ = E^j-EA' 

3=1 i=i 

B = (27r)(5-*'-^)/2^"/^^-i/2 f[(A,)-"^+°-^n(5,- 

j=i 1=1 



,6,-0.5 



In our case, p = 1, q = 1, ai = 1 — 7/2, Ai = 7/2, 61 = and Bi = 1. Therefore, for z large 



Hi:', 



(1-7/2,7/2) 

(0,1) 



V(2-7)^ 



7 



(7-l)/(2-7) 



,(7-l)/(2-7) 



exp 



2 _ 7/(2-7) 



2/(2-7) 



(3.36) 



Substituting this in the expression for f{t) [cf. Eq. ( p.22[ )], we obtain for small t 

d 



fit) 



where 



t(4-7)/(4-27) 



07 



exp 



t7/(2-7) 



(3.37) 



07 



(7-l)/(2-7) 



4ir7r(2 -7) V2/^/ 

2-7^7y/(^-^V « 

7 



For < 7 < 2, both r and c? are greater than zero. Therefore, f{t) is exponentially decaying 
for small t. Note that for 7 = 1 (ordinary Brownian motion), the above expression reduces 
to 



fit) 



exp[-a'^/4:Kt]. 



(3.38) 
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In this case, the asymptotic expression turns out to be the exact expression. 

We can determine the t value where f{t) attains its maximum value from the above 
expression. We obtain 

W = -r^ . (3.39) 



For ordinary Brownian motion (7 = 1), 



7, 



2d , ^ 



The theoretical prediction for the full FPT density function given in Eq. ( |3.22|) is verified 
by numerically simulating the underlying CTRW process characterized by the probability 
density function (j){y,u) [cf. Eq. ( p.5| )]. For the sake of numerical efficiency, we replace the 



waiting time distribution (a — 1)/t(1 + u/r)" in u) by the Pareto density function ||T6 
(which is equal to zero for u < t and (a — l)r"^^/u° for u > r). The parameter values are 
chosen as follows: 7 = 0.5, a = 1.0 and K = 0.1. In Figure 2, the FPT density function 



obtained theoretically from Eq. (|3.22|) is compared with the FPT density function obtained 



numerically using 10 million realizations of the underlying stochastic process. If we take a 
large value for r (=0.01), the agreement is not that good (see Figure 2a) since we are not 
yet close to the generalized diffusion limit. If, however, we take r to be 10"*^, the numerical 
simulation is in excellent agreement with the theoretical prediction (see Figure 2b). The 
agreement gets better as r and a become smaller, approaching the generalized diffusion 
limit. 

The above statement can be quantified as follows. Consider the Kullback-Leibler in- 
formation criterion ||2^j3^ which gives a quantitative measure of how "far apart" a given 



approximate density function fi{t) is from the exact density function f{t). The Kullback- 
Leibler information criterion is defined as 

KL{f, h) = r dtfit) log ^ (3.41) 

This is always greater than or equal to zero and is equal to zero if and only if /i(t) agrees 
exactly with f{t). The deviation away from zero quantifies the disagreement between the 
two density functions. 

In our case, we take fi{t) to be the approximate numerically simulated FPT density 
function for various values of r (the value of a is correspondingly varied to keep K constant). 
The plot of the Kullback-Leibler information criterion for different values of 1/r is given in 
Figure 3. As r decreases (i.e., as the generalized diffusion limit is approached), the Kullback- 
Leibler information criterion also decreases. 



B. Nonzero Drift Case 

In this section, we consider the first passage time problem for Levy type anomalous 
diffusion with drift. An explicit solution for the FPT density function is not possible in this 
case. We obtain the Laplace transform of the FPT density function (the so-called "moment 
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generating function" [0). Using this, we obtain the mean and variance of the first passage 
time distribution. For regular Brownian motion (7 = 1), the inverse Laplace transform of 
the moment generating function can be explicitly carried out to give the standard results. 
Taking the Laplace transform of Eq. (|2.33|) we get 



W{x,Q) K /i 9 

s ox^ s ox 



(3.42) 



where q{x,s) is the Laplace transform of W{x,t). Here we have again applied the result 
that the Laplace transform of oD^^W{x,t) is q{x,s)/s'^. The above equation can be 



rewritten as 



where 



d"^ d 

— g(x, s) + A—q{x, s) + Bq{x, s] 



K 



-6{x — a) 



A 



7-1 



K 



Since s,K > 0, we have 



2 27-2 7 



(3.43) 



(3.44) 



(3.45) 



Therefore two independent solutions of the homogeneous equation corresponding to Eq. 
( p.43| ) are given by |3^ 



gi(x, s) = exp[a;(A - A)/2]; q2{x, s) = exp[a;(-A - A)/2]. 



(3.46) 



Consequently, the general solution of Eq. ( |3.43| ) satisfying all the boundary and initial 
conditions [cf. Eq. ( p.2|) ] is given by 



q{x,s] 



.7-1 



KX 



_g -Aix-a)/2l^-\\x-a\/2 _ ^-\(x+a)/2l 



(3.47) 



To obtain the Laplace transform of the FPT density function, we take the Laplace 
transform of Eq. (|0|) to get 



F{s) = -s / dxq{x,s)+ / dxW{x,0). 
Jo Jo 



(3.48) 



Here we have used the fact that Laplace transform of dW{x, t)/dt is given by ||T3| sg(x, s) — 
W{x,^). Since W{x,^) = 5{x - a) [cf. Eq. (pj)], we obtain 



/•OO 

F{s) = 1 — s dx q{x, s) 
Jo 

Substituting for q{x, s) from Eq. ( p.47|) , we get 



(3.49) 



F{s) = 1 - 



KX 

s7 roD 



+ 



KX Jo 
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The integrals can be easily evaluated to finally give [upon using Eq. ( p.44| )] 



F{s) = exp 



a/is 



7-1 



2K 



— a\ 



/ ^2^27-2 



(3.50) 



For < 7 < 1, F{s) is not a completely monotone function [Q. That is, it does not 
satisfy the following conditions 



flF"'( s"! 

-1)"^^ > for s > 0; F(0) = 1. 



(3.51) 



Hence F{s) is not a Laplace transform of a probability density function. The physical reason 
for this is not yet understood. For 1 < 7 < 2, we have not been able to prove rigorously that 
-^7,/i(s) is a completely monotone function. However, we have calculated the first hundred 
derivatives of F^^^{s) using the symbolic manipulation program Mathematica and find that 
all of them satisfy Eq. ( |3.51| ). Hence we conjecture that F{s) is the Laplace transform 
of a probability density function for 1 < 7 < 2. Henceforth, we restrict ourselves to this 
parameter range. 

First consider the case where the drift is towards the barrier. This implies that /i < 
since in our formulation the diffusive process starts at x = a > and the barrier is at x = 0. 
In this case, the FPT density function can be written as [cf. Eq. (|3.5CI| )] 



F{s) 



where 



|/i|S 



7-1 



7-1 



2K 2K 
The mean first passage time is given by 

dF{s] 



1 + 



{T)- 

From Eqs. ( ^.52| ) and ( |3.53| ), we have 

dF{s) _ 
ds 



ds 



s=Q- 



(3.52) 



(3.53) 



(3.54) 



|/x|(7-l)g^-^a _ \ ^ 4ir.2-7 \ _ (2-7Hl + ^)-^/^ 



2K 



(3.55) 



We need to find the limiting value of the above expression as s 
As s — > 0, we can expand the square root in Eq. ( p.53| ) to give 

s Ks^-^ 
g{s) = -^ + ^ 



0. First consider e"^*^*-*. 



(3.56) 



Hence g{s) ^ as s — 0. Consequently, e'^^'-'^-* — 1 as s — 0. Performing similar expansions 
for the other terms in Eq. ( p.55|) , we finally obtain [cf. Eq. (|3.54| )] 
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Note that the mean first passage time is independent of 7. 
The variance is obtained as follows: 

Var(T) = (T^) - {Tf. 

Therefore we need to evaluate (T^). This is given by 

(PF{s) 



in 



s=0- 



Now 



(PF{s) 



d^g{s) ( dg{s) 
a — — h a- 



ds'^ ' \ ds 

Consider the first term. The second derivative of g{s) is given by 



ds'^ 



2K 

(7-l)(2-7) 



1 + 



'1/2 



-3/2 



Expanding all terms we get 



ds^ 



|/i|(7 - 1)(7 - 2)s^-3 ^ . 3 ■ • ■ (2n - 3) (AKs"^ 



2K 



E 

n=2 



2"n! 



(7 - 1)(2 - 7) ^ ■ 3 ■ ■ ■ (2n - 1) /4irs 



E 

n=l 



2"n! 



- o/^2„l-7 00 



2ir(2-7)^s 



E 

n=0 



;-!)"! • 3 ■ ■ ■ (2n + 1) /4irs2-T\" 



13 ^ 2"n! 

After considerable manipulation, this can be rewritten as 



d^g{s) _2K{2-^)s'^^ 



ds^ 



n=0 L 



n(2-7) + (3-7) 



(n + 2) 



-1)"1 ■3---(2n + 1) (AKs 



2-7^ 



2"n! 



Thus the first term in Eq. ( p.60|) is given by 



ds^ 



ag{s) 



E 



n=0 L 



n(2-7) + (3-7)' 



{n + 2) 



-l)"l-3---(2n + l) fAKs'^-^' 



2"ra! 
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For 7 = 1, the case of ordinary Brownian motion, we obtain from the above equation 



A(£) 



s=0 



(3.62) 



Here we have also used the fact that K = for 7 = 1. On the other hand, for 1 < 7 < 2, 
as s — i> the prefactor multiplying the sum in Eq. ( p.61| ) diverges whereas the sum itself is 
finite and bounded away from zero. Consequently, for 1 < 7 < 2 the first term in Eq. (p. 60 ) 
diverges as s — 0. 

Next, we consider the second term in Eq. ( p.60|) . We already know the limiting behaviour 
of this term for l<7<2ass— >0 from our earlier analysis for mean first passage time, 
namely, 

2 ^2 

(3.63) 



'dg{s) 



ds 



s=0 



Substituting the above results in Eq. (|3.59|) , for 7 = 1 we obtain 



1/^ 



a 

3 ^ Ji?' 



Therefore [cf. Eq. ( ^381) ] 



Var(T) = (T^a/|/i| 



(3.64) 



(3.65) 



for 7 = 1. For 1 < 7 < 2, the first term in Eq. ( p.60D diverges as s ^ whereas the second 
term is finite. Therefore (T^) diverges. Hence the variance also diverges. 

Next, consider the case where the drift is away from the barrier. This implies that /i > 
in our formulation. Now the FPT density function can be written as in Eq. ( p.52|) where 



7-1 



7-1 



2K 



2K 



'1 + 



(3.66) 



Performing the same analysis as above, it is easily seen that the mean and variance diverge 
for 1 < 7 < 2. 

For 7 = 1, the Laplace transform of the first passage time density function reduces to 
[cf. Eq. (pD] 



F(s) 



exp 



Now the inverse Laplace transform of exp(— ay^), a > is |T3 



a 



exp(-a^/4t). 



(3.67) 



(3.68) 



Using this result, we can easily perform the inverse Laplace transform of Eq. ( p.67| ) to obtain 



fit) 



exp 



(a + /it)" 
AKt 



a > 0, if: > 0. 



(3.69) 



We comment that, if the starting point of the diffusion is chosen at a;(0) = 0, a negative 
/i will be used in the above equation and we recover the expected inverse Gaussian density 
||T9| for the FPT density function of Brownian motion with drift . 
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C. Absorbing Barriers at x = 



—b and x = a 



In this section, we consider the first passage time problem for Levy type anomalous 
diffusion with zero drift and absorbing barriers placed at x = —b and x = a. We obtain an 
explicit solution for the first passage time density. 

The FFPE to be solved is given in Eq. ( |3.4| ). The boundary and initial conditions 
become: 



W{-b, t) = W{a, t) = 0, W{x, 0) = 6{x). 



(3.70) 



We solve the FFPE using the method of separation of variables ||2T|]. Let W{x, t) = X(x)T{t). 
Substituting in Eq. ( |3.4|) we obtain 

(3.71) 



X{x)T{t)-X{x) = oD;''T{t)KX"{x), 



where X"{x) denotes the second derivative of X{x) with respect to x. Separating out the 
variables and introducing the separation constant A we get 



KX"{x) = XX (x), 



and 



T{t) - 1 = A oD^^nt). 
The solution of Eq. ( |3.72| ) with the given boundary conditions is given by 

n7r(6 + x) 



XJx) = A„,sm 



{a + b) 



with 



A. 



f , m"^^' n = l,2,... . 
{a + by 



(3.72) 



(3.73) 



(3.74) 



(3.75) 



To solve Eq. (|3.73|) , we take its Laplace transform to obtain (introducing the subscript n 
coming from A„) 



Tn{s) 



1 A. 



'-TJs). 



(3.76) 



Here we have used the fact that the Laplace transform of qD^ '^T{t) is T[s)/s"'. Solving for 
TJs) we get 



TJs) 



(3.77) 



Taking the inverse Laplace transform [22| we finally obtain 



Tn{t) = E-y 



(a + by 



(3.78) 
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where E^{z) is the Mittag-Leffler function introduced earher [cf. Eq. (|3.8|) ]. 
Combining the solutions for space and time parts, we get 



sm 



n=l 



n7r{b + x) 
(a + b) 



(a + h) 



(3.79) 



The coefficients An are determined by imposing the initial condition W{x,Q) = 6{x). This 
gives us 



a + b 



sm 



Hence 



W{x,t) = J2 rsm 



nnb 
(a + b) 



sm 



mrb 
a + b) 

n7r{b + X 
(a + b) 



T'la + b 

From Eq. ( |3.3| ), the FPT density function is given by 

d 



(a + by 



-Kt^ 



(3.80) 



(3.81) 



fit) 



dt J-b 



dx W{x,t). 



(3.82) 



Substituting for W{x,t) into this equation, we get 



m 



sm 



(2ra + l)7r6 



dt [ (2n + l)7r 

Here we have used the fact that 

mi{b + x) 
{a + b) 



E^ 



f 

J-b 



sm 



[a + b) 

2(a + b) 



{2n + 1)V 
{a + by^ 



-Kt^ 



(3.83) 



if n is odd, 



riTC 

if n is even. 



To evaluate the derivative in Eq. ( 3.83|) , we write the Mittag-Leffler function in terms of 
the H-function using Eq. ( |3.11| ): 



fit) 



4 °° 

E 



sm 



(2n+ l)7r6' 
ia + b) 



d ,1 f{2n + 1)^71^ 



dt 



(a + 6)^ 



-Kt^ 



(0,1) 

(0,1), (0,7) 



(3.84) 



But [0 



d ^,,,fi2n + l)\' 



dt 



{a + by 



■Kt^ 



(0,1) 

(0,1), (0,7) 



1 / (2n + l)V 

t^^H i<^ + by 

1 / (2n + l)V 



(0,7) (0,1) 
(0,1), (0,7) (1,7) 

(0,1) \ 

(0,1), (1,7) r 



In the last step, we have used properties 1 and 2 of H-functions [cf. Appendix]. From 
the series expansion of the H-function given in Eq. (|3.10| ), it can be shown after some 
manipulation that 
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(0,1) 

(0,1), (1,7) 



where E^^piz) is the generahzed Mittag-Leffler function |^| 



Using this in the expression for f{t) we get 



> 0. 



/(^) = T—^E (2^ + 1) sin 

n=0 



(a + 6)2 



f{t) 



^ (2n + 1) sin 

n=0 



(a + 6)2 



"(2n + l)7r6' 




(a + 6) 


, we obtain 




"(2n+ l)7r6" 


exp 


(a + 6) 



{2n + l)2 7r2 
(a + 6)^ 



(2r2+ l)27 r2 
(a + 6)2 



(3.85) 



(3.86) 



(3.87) 



(3.88) 



IV. CONCLUSIONS 

We studied the first passage time (FPT) problem for Levy type anomalous diffusion. 
We obtained the FPT distribution using the recently formulated fractional Fokker-Planck 
equation. We derived an explicit expression for the FPT distribution in terms of Fox or H- 
functions when the diffusion has zero drift. The theoretical result was verified by numerically 
simulating the underlying continuous time random walk. When the drift is nonzero, we 
obtained an analytic expression for the Laplace transform of the FPT distribution. This 
was used to calculate the mean and variance of the FPT distribution. Finally, for the case of 
two absorbing barriers at finite distances from the origin we expressed the FPT distribution 
in terms of a power series. In all of the above situations, the known results for ordinary 
diffusion (Brownian motion) were obtained as special cases of our more general results. 
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APPENDIX: PROPERTIES OF H-FUNCTIONS 



The H-function has the following remarkable properties |]TD[ which we will use later. 
Property 1. The H-function is symmetric in the pairs (ai, Ai), . . . , (a„, A„), likewise 



, yap, Ap)] in {pi, Bi), . . . , (6^, B^^) and in {bm+i, ^m+i), • • • , Bq) 



Property 2. Provided n > 1 and q > m 



(ai,Ai), (02,^2), 
(&i,5i), 

(02,^2), ■ ■ - , (ap,v4p) \ 
(61, (6,-1, )• 



Property 3. Provided m > 2 and p > n 



Trm-l,n 
-"p-l,g-l I ^ 



(ai,Ai), (ap_i,y4p_i), 
(61, Si), (62,52), (6„5,) 

(ai,Ai), ■ ■ ■ , (ap_i, Ap_i) \ 
(62,52), (6„5,) )• 



Property 4- 



HI 



{bj,Bj)j=i^ 



TTn,m 



Property 5. For A; > 0, 



(flj, Aj)j=l^,,,^p 

(6j-, 5j)j=l,...,q 



Property 6. 



(Oj, Aj)j=l,...,p 



{h,B 



h: 



j, ^j)j = l,...,q 



p,g 



(1 - 6j,5j)j = l„„,g 

(1 — aj, Aj)j=i^ p 



(oj, A;Aj)j=i^...^p 
(6j, /c5j)j=i,...,g 



(Oj + pAj, Ajr)j = l,...,p 
{bj + pBj,Bj)j = i^,„^q 



(Al) 



(A2) 



(A3) 



(A4) 



(A5) 
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APPENDIX: FIGURE LEGENDS 



Figure 1: A log- log plot of the mean squared displacement obtained by numerically simu- 
lating the underlying CTRW process for a Levy type anomalous diffusion with 7 = 0.5. 

Figure 2: a) Comparison of the theoretical FPT distribution (solid line) with the distribu- 
tion (dashed line) obtained by numerically simulating the underlying CTRW process 
for a Levy type anomalous diffusion (with 7 = 0.5) that is not close to the generalized 
diffusion hmit (r = 10~^). b) Same comparison as above but closer to the generahzed 
diffusion limit (r = 10~^). 

Figure 3: Plot of KuUback-Leibler information criterion against l/r. 
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